c=========================================================================
C  The following routine solves the differential equation below backwards
C  from S(JSTOP) = S3 and S(JSTOP-1) = S2
C        2
C       D S   [ LN(LN+1)       V(X)          ]
C       --- - [ --------  + --------- - ECMN ] S  = 0
C         2   [   X*X           X            ]
C       DX      
C
C  Note that V(X)/X = UCENTR(X)
      subroutine nmrvb8(ln,ecmn,ucentr,cntfug,gridx,nx,
     >   jdouble,njdouble,s2,s3,chi,jstart,jstop)
c
      implicit real*8 (a-h, o-z)
c
      dimension ucentr(nx),cntfug(nx),gridx(nx),jdouble(njdouble),
     >   chi(nx)

      j = jdouble(njdouble-1)
      dx =  gridx(j+1)-gridx(j)
      h2 =  dx * dx
      h2d =  h2/12d0
      xl =  gridx(jstop)
      xlm1 =  gridx(jstop-1)
      wnn = sqrt(abs(ecmn))

      f3 =  ucentr(jstop)+cntfug(jstop)-ecmn
      f2 =  ucentr(jstop-1)+cntfug(jstop-1)-ecmn
      t3 =  (1d0 -h2d * f3) * s3
      t2 =  (1d0 -h2d * f2) * s2

      chi(jstop) =  s3
      chi(jstop-1) =  s2

      istart = njdouble-1
      do while (istart.gt.1.and.jdouble(istart).gt.jstop)
         istart = istart-1
      end do
      istop = istart
      istart = istart+1
      do while (istop.gt.1.and.jdouble(istop).gt.jstart)
         istop = istop-1
      end do
      istop = istop+1
C  JDOUBLE(ISTOP) points to the first doubling of DX that happens after JSTOP
C  JDOUBLE(ISTART) points to the last doubling of DX that happens before JSTART
      do i = istart,istop,-1
         j1 = min(jstop,jdouble(i))-2
         j2 = max(jstart,jdouble(i-1))
         do j = j1,j2,-1
            f1 =  ucentr(j)+cntfug(j)-ecmn
            t1 =  2d0 * t2 +h2 * f2 * s2 -t3
            s1 =  t1/(1.0d0- h2d * f1)
            t3 = t2
            t2 = t1
            f3 = f2
            f2 = f1
            s3 = s2
            s2 = s1
            chi(j) =  s1
         end do
         if (j2.eq.jstart) return
         j = j2-1
         dx =  dx/2d0
         h2 =  dx * dx
         h2d =  h2/12d0
         f1 =  ucentr(j)+cntfug(j)-ecmn
         s1 =  s2 * (36d0 +33d0 * h2 * f2) +s3 * (-12d0+5d0 * h2 * f3)
         s1 =  s1/(24d0 +2d0 * h2 * f1)
         t2 =  s2 * (1d0 -h2d * f2)
         t1 =  s1 * (1d0 -h2d * f1)
         t3 = t2
         t2 = t1
         f3 = f2
         f2 = f1
         s3 = s2
         s2 = s1
         chi(j) = s1
      end do
      
      return
      end
